Controle de processos βAR(1)

# Temas para os gráficos

theme.base <- theme_minimal(base_size = 11) +
  theme(
    axis.text = element_text(size = 8),
    plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 9),
    plot.caption = element_text(size = 8),
    axis.title = element_text(size = 8),
    legend.title = element_text(size = 8),
    panel.grid.major = element_line(colour = "grey90", linewidth = 0.5),
    panel.grid.minor = element_line(
      colour = adjustcolor("grey90", alpha.f = 0.4),
      linewidth = 0.5
    ),
    panel.border = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    axis.line.x = element_line(colour = "grey"),
    axis.line.y = element_line(colour = "grey"),
  )

theme.no_legend <- theme(legend.position = "none")

theme.no_grid <-  theme(panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank())

theme.no_axis <- theme(axis.line.x = element_blank(), axis.line.y = element_blank())

theme.no_axis_title <- theme(axis.title.x = element_blank(), axis.title.y = element_blank())

# Theme for timeseries
theme.ts <- theme.base +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank()
  ) +
  theme.no_legend

plotly.base <- function(p) {
  p %>%
    layout(margin = list(b = 60, t = 80)) %>%
    config(mathjax = 'cdn')
}

# Função para salvar e carregar resultados
cache_dados <- function(chave, funcao_geradora) {
  c <- paste("cache", chave, sep = "/")
  c <- paste(c, "rds", sep = ".")
  cache <- file.exists(c)
  
  if (!cache) {
    resultado <- funcao_geradora()
    saveRDS(resultado, c)
  } else {
    resultado <- readRDS(c)
  }
  
  resultado
}

Introdução

O βAR(1) é um modelo de séries temporais autorregressiva de ordem 1 com distribuição Beta – onde a variável de interesse está restrita no intervalo (0, 1) –, este modelo foi proposto por Rocha e Cribari-Neto:

• Beta autoregressive moving average models. 2009. DOI: https://doi.org/10.1007/s11749-008-0112-z
• Erratum to: Beta autoregressive moving average models. 2017. DOI: https://doi.org/10.1007/s11749-017-0528-4

Sejam a série temporal \(Y_{t} = \left\{y_{1}, y_{2}, \ldots, y_{t} | y_{n} \in (0, 1)\right\}\), o conjunto \(l\)-dimensional de covariáveis \(\mathbf{X_{t}}\), e o vetor de parâmetros \(\beta = \left\{\beta_{1}, \beta_{2}, \ldots, \beta_{l}\right\}\), temos o modelo βAR(1) é definido como:

\[ g(\mu_t) = \alpha + \mathbf{X_t'}\beta + \phi \left(g(Y_{t-1}) - \mathbf{X_{t-1}\beta}\right) \]

Onde \(g: (0,1) \rightarrow \mathbb{R}\) é a função de ligação, \(\alpha\) é o intercepto, \(\phi\) é o parâmetro de autorregressão, e \(\mu_t\) é a média condicional da distribuição Beta. A função de ligação mais usual é a função logit.

Neste trabalho, objetiva-se aplicar métodos de controle de processos βAR(1) para detectar mudanças no processo. Para isso, utilizamos que, quando modelada corretamente, os resíduos de uma série temporal, como a βAR(1), são independentes e normalmente distribuídos com média zero e variância constante.

Desta forma, quando o processo sofre uma mudança, espera-se que os resíduos não sejam mais independentes e que a média e a variância dos resíduos sejam diferentes. Assim, podemos utilizar métodos de controle de processos para detectar essas mudanças.

Modelo

Olhando sob a perspectiva de teste de hipóteses, temos:

• \(H_0\): amostra inicial de tamanho 100, com \(\Phi_0 = 0.2\). Indicando a Fase I da aplicação do controle de processos.
• \(H_1\): amostra subsequente de tamanho 200, com \(\Phi_1 = 0.2, 0.3, \ldots, 0.6\). Indicando a Fase II da aplicação do controle de processos.

Além disso, utilizaremos um nível de significância de 0.05 para os testes de hipóteses, o que nos dá um quantil de 1.96 para a distribuição normal padrão.

A simulação dos dados será feita a partir da biblioteca BTSR: Bounded Time Series Regression.

nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2
alpha_teste <- 0.05
quantil_teste <- qnorm(1 - alpha_teste / 2)

coeficientes <- function(phi) {
  # βARMA: o modelo de Rocha e Cribari-Neto (2009, 2017)
  #        é obtido definindo `coefs$d = 0`
  #        e `d = FALSE` e `error.scale = 1` (escala preditiva)
  #
  # ν (nu): parâmetro de precisão, quanto maior o seu valor,
  #         menor a variância condicional (para μ_t fixo).
  #         O pacote BTSR define como padrão ν = 20.
  list(
    alpha = 0,
    nu = 20,
    phi = phi,
    d = 0,
    nu = 20
  )
}

barma.sim <- function(n,
                      phi,
                      seed = NULL,
                      y.start = NULL) {
  BARFIMA.sim(
    n = n,
    coefs = coeficientes(phi),
    y.start = y.start,
    error.scale = 1,
    complete = F,
    seed = ifelse(is.null(seed), sample(1:1000, 1), seed)
  )
}

barma.phi_estimado <- function(yt,
                               alpha = 0,
                               nu = 20,
                               phi = 0.1) {
  BARFIMA.fit(
    yt = yt,
    start = list(alpha = alpha, nu = nu, phi = phi),
    p = 1,
    # Odem do polinômio AR
    d = FALSE,
    error.scale = 1,
    report = F
  )$coefficients["phi"][[1]]
}

barma.residuos <- function(yt, phi_estimado) {
  BARFIMA.extract(
    yt = yt,
    coefs = coeficientes(phi_estimado),
    llk = F,
    info = F,
    error.scale = 1
  )$error
}

# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func,
                                  debug = F,
                                  max = 5,
                                  ...) {
  FINALMENTE <- F
  contador <- 0
  
  while (!FINALMENTE) {
    contador <- contador + 1
    if (debug)
      print(paste("Tentativa:", contador))
    
    resultado <- tryCatch({
      func(...)
    }, error = function(err) {
      if (contador >= max) {
        stop("Deu ruim, número máximo de tentativas atingido")
      }
      
      if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
        # Ignora erro de extração de resíduos
        return(NULL)
      }
      
      stop(err)
    })
    
    if (!is.null(resultado)) {
      FINALMENTE <- T
    }
  }
  
  return(resultado)
}

Monte Carlo

# Função auxiliar para gerar os dados
# Exemplo 1:
# ```r
# teste1 <- gerador_monte_carlo(
#   parametros = list(
#     lambda = seq(0.1, 0.4, by = 0.1)
#   ),
#   numero_de_execucoes = 2
# )
#
# teste1
# ```
#
# Exemplo 2:
# ```r
# teste2 <- gerador_monte_carlo(
#   parametros = expand.grid(
#     desvio_detectavel = seq(0.6, 1.8, by = 0.2),
#     intervalo_de_decisao = seq(4, 6, by = 1)
#   )
#   numero_de_execucoes = 2
# )
#
# teste2
# ```
#
# Total de execuções é:
# `execucoes = numero_de_execucoes * length(novos_phis) * nrow(parametros)`
#
gerador_monte_carlo <- function(parametros = NULL,
                                numero_de_execucoes = 200) {
  # Verifica se `parametros` é um data frame, caso contrário converte para um data frame
  if (!is.null(parametros) && !is.data.frame(parametros)) {
    parametros <- as.data.frame(parametros)
  }
  
  novos_phis <- seq(0.2, 0.6, by = 0.1)
  
  # Expande a grade de parâmetros para cobrir todas as combinações
  result <- expand.grid(k = 1:numero_de_execucoes, h1_phi = novos_phis)
  
  if (!is.null(parametros)) {
    result <- merge(result, parametros, all = TRUE)
  }
  
  result %>%
    mutate(id = row_number()) %>%
    rowwise() %>%
    mutate(
      # H0
      ## Gera amostra de controle
      h0_amostras = list(barma.sim(
        n = nH0, phi = phi_parametro, seed = id
      )),
      ## Estima o valor de phi da amostra de controle
      h0_phi = barma.phi_estimado(h0_amostras),
      ## Calcula os resíduos da amostra de controle
      h0_residuos = list(barma.residuos(h0_amostras, h0_phi)),
      # H1
      ## Gera a amostra subsequente
      h1_controle = list(
        barma.sim(
          n = nH1,
          phi = h1_phi,
          # última observação da amostra de controle
          y.start = h0_amostras[nH0],
          seed = id + 1337E4
        )
      ),
      ## Calcula os resíduos da amostra subsequente
      h1_residuos = list(so_quero_que_funcione(\() barma.residuos(
        h1_controle, h0_phi
      )))
    )
}

Validação

Vamos verificar nossa implementação do gerador de Monte Carlo.

teste_montecarlo <- cache_dados("teste-montecarlo", function() {
  gerador_monte_carlo(parametros = list(nH0 = c(
    rep(25, 20), rep(50, 20), rep(100, 20), rep(200, 20)
  )),
  numero_de_execucoes = 1) %>%
    select(nH0, h0_phi, h1_phi) %>%
    mutate(# Calcula a diferença entre os valores de phi
      # `phi_parametro`: valor de phi definido para a amostra de controle
      # `h0_phi`: valor de phi estimado para a amostra de controle. Espera-se que seja igual a `phi_parametro`
      diferenca = h0_phi - phi_parametro)
})

teste_montecarlo %>%
  group_by(nH0) %>%
  summarise(
    simulações = n(),
    mean = mean(diferenca),
    var = var(diferenca),
    min = min(diferenca),
    max = max(diferenca),
    .groups = "drop"
  )

ggplotly(
  teste_montecarlo %>%
    ggplot(aes(x = factor(nH0), y = diferenca)) +
    geom_boxplot() +
    geom_hline(
      yintercept = 0,
      color = "red",
      linetype = "dashed"
    ) +
    annotate(
      geom = "text",
      x = 0.55,
      y = 0 + 0.01,
      label = "0.00",
      color = "red"
    ) +
    labs(
      x = "Número de execuções",
      y = "Diferença entre Φ e Φ₀",
      title = paste0("Diferença entre Φ e Φ₀")
    ) +
    theme.base
) %>%
  plotly.base

EWMA

O EWMA (Exponential Weighted Moving Average) é um método de controle de processos que utiliza uma média móvel ponderada exponencialmente para detectar mudanças no processo. Sendo definido por:

\[ z_i = \lambda y_i + (1 - \lambda) z_{i-1} \]

Onde, \(z_i\) é a estatística de controle no instante \(i\), \(y_i\) é a observação no instante \(i\), \(\lambda\) é o fator de suavização e \(z_{i-1}\) é a estatística de controle no instante anterior.

O fator \(\lambda\) é uma constante definida no intervalo \((0, 1]\) e seu valor inicial (para \(i = 1\)) é definido como a média do processo, tal que \(z_0 = \mu_0\).

A estatística de controle, \(z_i\), é comparada com os limites de controle, \(\text{LCS}\) e \(\text{LCI}\), sendo é definido como:

\[ \text{LC} = \mu_0 \pm L \sigma \sqrt{\frac{\lambda}{2 - \lambda} \left[1 - \left(1 - \lambda\right)^{2i}\right]} \]

Aqui, utilizamos o pacote qcc para implementar o EWMA.

ewma_qcc <- function(amostra_inicial,
                     lambda,
                     amostra_subsequente,
                     nsigmas = 1.96) {
  ew <- ewma(
    amostra_inicial,
    lambda = lambda,
    nsigmas = nsigmas,
    plot = F,
    newdata = amostra_subsequente
  )
  
  registros <- (nH0 + 1):(nH0 + nH1)
  
  ewma <- as.numeric(ew$y[registros])
  LI <- ew$limits[, 1][registros]
  LS <- ew$limits[, 2][registros]
  fora_de_controle <- ewma < LI | ewma > LS
  total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
  fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
  list(
    ewma = ewma,
    LI = LI,
    LS = LS,
    fora_de_controle = fora_de_controle,
    total_fora_de_controle = total_fora_de_controle,
    fracao_fora_de_controle = fracao_fora_de_controle
  )
}

Monte Carlo para α = 0.05

Vamos analisar buscar o melhor valor de \(\lambda\) para o EWMA, mantendo um nível de significância constante, \(\alpha = 0.05\).

ewma_monte_carlo <- cache_dados("simulacao-ewma", function() {
  gerador_monte_carlo(parametros = list(lambda = seq(0.1, 0.9, by = 0.1))) %>%
    mutate(
      qcc = list(ewma_qcc(h0_residuos, lambda, h1_residuos)),
      ewma = list(qcc$ewma),
      LI = list(qcc$LI),
      LS = list(qcc$LS),
      fora_de_controle = list(qcc$fora_de_controle),
      total_fora_de_controle = qcc$total_fora_de_controle,
      fracao_fora_de_controle = qcc$fracao_fora_de_controle
    ) %>%
    select(-qcc)
})

ewma_monte_carlo$lambda <- as.factor(ewma_monte_carlo$lambda)

Cartas para λ = 0.2

Vamos analisar o EWMA para \(\lambda = 0.2\). Aqui, vamos considerar apenas a primeira execução.

ggplotly(
  ewma_monte_carlo %>%
    filter(k == 1 & lambda == 0.2) %>% # Apenas a primeira execução
    select(h1_phi, ewma, LI, LS, fora_de_controle, lambda) %>%
    rename(`Φ₁` = h1_phi) %>%
    mutate(observacao = list(1:nH1)) %>%
    unnest(cols = c(
      `Φ₁`, observacao, ewma, LI, LS, fora_de_controle
    )) %>%
    ggplot() +
    geom_line(aes(
      x = observacao, y = LI, color = "Limite Inferior"
    ), linetype = "dashed") +
    geom_line(aes(
      x = observacao, y = LS, color = "Limite Superior"
    ), linetype = "dashed") +
    geom_line(aes(
      x = observacao, y = ewma, color = "EWMA"
    )) +
    geom_point(
      data = . %>% filter(fora_de_controle),
      aes(x = observacao, y = ewma, color = "Fora de controle"),
      size = 1,
      shape = 4
    ) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
        "Limite Superior" = adjustcolor("red", alpha.f = 0.5),
        "EWMA" = adjustcolor("black", alpha.f = 0.8),
        "Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.8)
      )
    ) +
    labs(title = "EWMA para resíduos βAR(1), com λ = 0.2", ) +
    theme.base +
    theme(legend.position = "bottom", legend.title = element_blank()) +
    facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
  plotly.base

Resumo

Resumo da fração de pontos fora de controle para diferentes valores de \(\lambda\) e \(\Phi_1\).

ewma_monte_carlo_resumo <- ewma_monte_carlo %>%
  group_by(lambda, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

datatable(
  ewma_monte_carlo_resumo %>%
    mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e λ",
  colnames = c("Valores de λ", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)

Gráficos

Pela análise gráfica das médias de FPFCs (Fração de Pontos Fora de Controle) notamos que \(\lambda = 0.7\) é o melhor valor para detectar mudanças no processo.

ggplotly(
  ewma_monte_carlo_resumo %>%
    ggplot(aes(
      x = h1_phi, y = mean, color = lambda
    )) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      color = "Valores de λ",
      title = "FPFCs por Φ e λ"
    ) +
    theme.base
) %>%
  plotly.base

E, a distribuição dos PFCs.

ggplotly(
  ewma_monte_carlo %>%
    ggplot(aes(
      x = factor(h1_phi),
      y = fracao_fora_de_controle,
      fill = lambda
    )) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      fill = "Valores de λ",
      title = "FPFCs por Φ e λ"
    ) +
    theme.base
) %>%
  plotly.base %>%
  layout(boxmode = 'group')

EWMA-AR

Segundo Montgomery (2009), o EWMA-AR é uma extensão do EWMA que utiliza um modelo autorregressivo para prever a próxima observação.

Assim, temos que \(\lambda \in (0, 1]\), sendo que, a previsão para a observação \(x_{t+1}\) é dada por \(\hat{x}_{t+1}(t)=z_{t} = \lambda x_t + (1 - \lambda) z_{t-1}\).

Otimizando λ

De acordo com Montgomery (2009), podemos encontrar um valor ótimo para \(\lambda\) através da minimização da soma dos quadrados dos resíduos.

E, ainda, temos que os erros de previsão são dados por \(e_t = x_t - \hat{x}_t(t-1)\) conforme a Eq. 10.16 (Montgomery, 2009).

Ou, de outra forma

\[ \hat{\lambda} = \min \left( \underset{\lambda\in(0,1]}{\arg\min}\,\text{Err}(\lambda) \right) \]

onde \(\hat{\lambda}\) é o \(\lambda\) ótimo por simulação, e \(\text{Err}(\lambda) = \sum_{t=1}^{n} e_t^2\).

Segundo, também Montgomery (2009), podemos estimar o valor de \(\sigma^2\) para o modelo EWMA-AR como \(\sigma^2 = \frac{\sum (\text{err}_i^2|_\lambda)}{n}\). Onde \(\text{err}|_\lambda\) são os resíduos do modelo EWMA-AR para o melhor valor de \(\lambda\) encontrado.

# Função para computar resíduos do EWMAAR
ewma_ar_residuos <- function(dados, ewma) {
  n <- length(dados)
  residuos <- numeric(n)
  residuos[1] <- dados[1]
  for (i in 2:n) {
    residuos[i] <- dados[i] - ewma[i - 1]
  }
  return(residuos)
}

# Função para computar o EWMA-AR
ewma_ar <- function(dados,
                    lambda,
                    x0 = NULL,
                    desvio = NULL) {
  desv <- ifelse(is.null(desvio), sqrt(sd(dados)), desvio)
  n <- length(dados)
  final <- ifelse(is.null(x0), n, n + 1)
  
  ewma_serie <- numeric(final)
  ewma_serie[1] <- ifelse(is.null(x0), dados[1], x0)
  for (i in 2:final) {
    ewma_serie[i] <- lambda * dados[i] + (1 - lambda) * ewma_serie[i - 1]
  }
  
  ewma <- ewma_serie[ifelse(is.null(x0), 1, 2):final]
  LI <- ewma - 1.65 * desv
  LS <- ewma + 1.65 * desv
  fora_de_controle <- dados < LI | dados > LS
  total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
  fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
  
  return(
    list(
      ewma = ewma,
      residuos = ewma_ar_residuos(dados, ewma_serie),
      LI = LI,
      LS = LS,
      fora_de_controle = fora_de_controle,
      total_fora_de_controle = total_fora_de_controle,
      fracao_fora_de_controle = fracao_fora_de_controle
    )
  )
}

lambda_otimo_optimize <- function(amostra_inicial) {
  lambda <- NA
  minimo <- Inf
  for (x in c(seq(0.01, 0.1, by = 0.01), seq(0.1, 0.9, by = 0.1))) {
    ew <- ewma_ar(amostra_inicial, lambda = x)
    
    erros <- ew$residuos
    soma_quadrados <- sum(erros ^ 2)
    if (soma_quadrados < minimo) {
      minimo <- soma_quadrados
      lambda <- x
    }
  }
  return(list(lambda = lambda, minimo = minimo))
}

ewmar_lambda_otimo <- cache_dados("ewmaar-lambda-otimo", \() {
  data.frame(id = 1:100) %>%
    rowwise() %>%
    mutate(
      h0_amostra = list(barma.sim(n = nH0, phi = phi_parametro)),
      h0_phi = list(barma.phi_estimado(h0_amostra)),
      h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
      ewma_ar = list(lambda_otimo_optimize(h0_residuos)),
      lambda = ewma_ar$lambda,
      minimo = ewma_ar$minimo
    ) %>%
    select(-ewma_ar)
})

ewmar_lambda_otimo %>%
  group_by() %>%
  summarise(
    mean_lambda = mean(lambda),
    mean_err = mean(minimo),
    var_err = var(minimo),
    min_err = min(minimo),
    max_err = max(minimo),
    .groups = "drop"
  )

Por simulação encontramos um \(\lambda\) ótimo de \(\sim 0.07\). Com erro médio de \(\sim 22\), o que nos dá um desvio de \(22 / 100 \simeq 0.2\).

Monte Carlo

ewmaar_monte_carlo <- cache_dados("simulacao-ewma-ar", function() {
  gerador_monte_carlo(parametros = list(lambda = seq(0.01, 0.1, by = 0.01))) %>%
    mutate(
      ewma_ar = list(ewma_ar(
        h1_controle, lambda, h0_amostras[nH0], sd(h0_amostras)
      )),
      ewma = list(ewma_ar$ewma),
      LI = list(ewma_ar$LI),
      LS = list(ewma_ar$LS),
      fora_de_controle = list(ewma_ar$fora_de_controle),
      total_fora_de_controle = ewma_ar$total_fora_de_controle,
      fracao_fora_de_controle = ewma_ar$fracao_fora_de_controle
    ) %>%
    select(-ewma_ar)
})

ewmaar_monte_carlo$lambda <- as.factor(ewmaar_monte_carlo$lambda)

Cartas para λ = 0.07

Vamos analisar o EWMA-AR para \(\lambda = 0.07\). O resultado encontrado anteriormente.

ggplotly(
  ewmaar_monte_carlo %>%
    filter(k == 1 & lambda == 0.07) %>% # Apenas a primeira execução
    select(h1_phi, h1_controle, LI, LS, fora_de_controle, lambda) %>%
    rename(`Φ₁` = h1_phi) %>%
    mutate(observacao = list(1:nH1)) %>%
    unnest(cols = c(
      `Φ₁`, observacao, h1_controle, LI, LS, fora_de_controle
    )) %>%
    ggplot() +
    geom_line(aes(
      x = observacao, y = LI, color = "Limite Inferior"
    ), linetype = "dashed") +
    geom_line(aes(
      x = observacao, y = LS, color = "Limite Superior"
    ), linetype = "dashed") +
    geom_line(aes(
      x = observacao, y = h1_controle, color = "Série"
    )) +
    geom_point(
      data = . %>% filter(fora_de_controle),
      aes(x = observacao, y = h1_controle, color = "Fora de controle"),
      size = 1,
      shape = 4
    ) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
        "Limite Superior" = adjustcolor("red", alpha.f = 0.5),
        "Série" = adjustcolor("black", alpha.f = 0.8),
        "Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.6)
      )
    ) +
    labs(title = "EWMA para resíduos βAR(1), com λ = 0.08", ) +
    theme.base +
    theme(legend.position = "bottom", legend.title = element_blank()) +
    facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
  plotly.base

Resumo

Resumo da fração de pontos fora de controle para diferentes valores de \(\lambda\) e \(\Phi_1\).

ewmaar_monte_carlo_resumo <- ewmaar_monte_carlo %>%
  group_by(lambda, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

datatable(
  ewmaar_monte_carlo_resumo %>%
    mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e λ",
  colnames = c("Valores de λ", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)

Gráficos

Pela análise gráfica das médias de FPFCs (Fração de Pontos Fora de Controle) notamos que o \(\lambda\) encontrado por otimização ultrapassa o limite de 5% de pontos fora de controle para \(\Phi_1 = 0.2\).

E, de uma forma geral, possui um resultado menos satisfatório que o EWMA sobre os resíduos.

ggplotly(
  ewmaar_monte_carlo_resumo %>%
    ggplot(aes(
      x = h1_phi, y = mean, color = lambda
    )) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      color = "Valores de λ",
      title = "FPFCs por Φ e λ"
    ) +
    theme.base
) %>%
  plotly.base

E, a distribuição dos PFCs.

ggplotly(
  ewmaar_monte_carlo %>%
    ggplot(aes(
      x = factor(h1_phi),
      y = fracao_fora_de_controle,
      fill = lambda
    )) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      fill = "Valores de λ",
      title = "FPFCs por Φ e λ"
    ) +
    theme.base
) %>%
  plotly.base %>%
  layout(boxmode = 'group')

CUSUM

O CUSUM (Cumulative Sum) é um método de controle de processos que utiliza a soma acumulada dos resíduos para detectar mudanças no processo.

Neste método são definidas duas constantes, \(H\) e \(K\), que representam o tamanho da mudança que se deseja detectar e a sensibilidade do método, respectivamente. As estatísticas de controle, sendo duas, são definidas como:

\[ \begin{matrix} \text{C}^+_i & = & \max\left[0, x_i - (\mu_0 + K) + C^+_{i-1}\right] \\ \text{C}^-_i & = & \max\left[0, (\mu_0 - K) - x_i + C^-_{i-1}\right] \\ \end{matrix} \]

onde \(x_i\) é a observação no instante \(i\), \(\mu_0\) é a média do processo, \(C^+_0 = C^-_0 = 0\) e \(i = 1, 2, \ldots, n\).

A constante \(K\), chamada de valor de referência, geralmente, segundo Montgomery (2009), é definida como a metade entre o \(\mu_0\) alvo e o valor fora de controle de média \(\mu_1\) que queremos detectar.

Já a constante \(H\), chamada de limite de decisão, é definida como o valor que a estatística de controle deve atingir para podermos detectar a mudança no processo. Para Montgomery (2009), um valor razoável para \(H\) é \(5\sigma\). Assim, se \(\text{C}^+_i \geq H\) ou \(\text{C}^-_i \leq -H\), então o processo está fora de controle.

cusum_qcc <- function(amostra_inicial_,
                      desvio_detectavel = 1,
                      intervalo_de_decisao = 5,
                      amostra_subsequente = NULL) {
  arguments <- list(
    data = amostra_inicial_,
    se.shift = desvio_detectavel,
    decision.interval = intervalo_de_decisao,
    plot = F
  )
  
  if (!is.null(amostra_subsequente)) {
    arguments$newdata <- amostra_subsequente
  }
  
  cu <- do.call(cusum, arguments)
  
  registros <- if (is.null(amostra_subsequente)) {
    seq(1, nH0)
  } else {
    seq(nH0 + 1, nH0 + nH1)
  }
  pos <- cu[["pos"]][registros]
  neg <- cu[["neg"]][registros]
  fora_de_controle_pos <- pos > intervalo_de_decisao
  fora_de_controle_neg <- neg < -intervalo_de_decisao
  total_fora_de_controle <- sum(ifelse(
    fora_de_controle_pos,
    fora_de_controle_pos,
    fora_de_controle_neg
  ))
  fracao_fora_de_controle <- total_fora_de_controle / length(registros)
  list(
    pos = pos,
    neg = neg,
    fora_de_controle_pos = fora_de_controle_pos,
    fora_de_controle_neg = fora_de_controle_neg,
    total_fora_de_controle = total_fora_de_controle,
    fracao_fora_de_controle = fracao_fora_de_controle
  )
}

Estimando o valor ótimo para K

Queremos que, sob \(H_1\), nas 100 amostras iniciais a probabilidade de um ponto fora de controle seja de, nominalmente, 5%. Portanto, vamos buscar um valor para K onde isso acontece.

desvio_detectavel_otimo_optimize <- function(amostra_inicial) {
  for (x in seq(0, 3, by = 0.1)) {
    fora <- cusum_qcc(amostra_inicial, desvio_detectavel = x)$total_fora_de_controle
    if (fora == 5) {
      return(x)
    }
  }
  return(NA)
}

cusum_k_otimo <- cache_dados("cusum-k-otimo", \() {
  data.frame(id = 1:1000) %>%
    rowwise() %>%
    mutate(
      h0_amostra = list(barma.sim(
        n = nH0, phi = phi_parametro, seed = id
      )),
      h0_phi = list(barma.phi_estimado(h0_amostra)),
      h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
      minimo = desvio_detectavel_otimo_optimize(h0_residuos)
    )
})

ggplotly(
  cusum_k_otimo %>%
    filter(!is.na(minimo)) %>%
    ggplot(aes(
      x = 0,
      y = minimo,
      group = 0,
      fill = "red"
    )) +
    geom_boxplot() +
    labs(x = element_blank(), y = "Valor ótimo de K", title = "Valor ótimo de K por simulação") +
    theme.base +
    theme.no_legend
) %>%
  plotly.base

valor_otimo_k <- mean(cusum_k_otimo$minimo, na.rm = TRUE)

kable(
  cusum_k_otimo %>%
    filter(!is.na(minimo)) %>%
    group_by() %>%
    summarise(
      "Média" = mean(minimo),
      "Desvio padrão" = sd(minimo),
      "Amostras" = n(),
    ),
  caption = "Valor ótimo para K",
  align = 'c',
  digits = 3
)

Valor ótimo para K

Média	Desvio padrão	Amostras
0.507	0.193	260

Monte Carlo

cumsum_k_monte_carlo <- cache_dados("simulacao-cusum-k", function() {
  gerador_monte_carlo(parametros = list(desvio_detectavel = seq(0.1, 1.0, by = 0.1))) %>%
    mutate(
      qcc = list(
        cusum_qcc(
          h0_residuos,
          desvio_detectavel = desvio_detectavel,
          amostra_subsequente = h1_residuos
        )
      ),
      pos = list(qcc$pos),
      neg = list(qcc$neg),
      fora_de_controle_pos = list(qcc$fora_de_controle_pos),
      fora_de_controle_neg = list(qcc$fora_de_controle_neg),
      total_fora_de_controle = qcc$total_fora_de_controle,
      fracao_fora_de_controle = qcc$fracao_fora_de_controle
    ) %>%
    select(-qcc)
})

Cartas para K ótimo

Vamos analisar o CUSUM para o valor ótimo de K encontrado.

cusum_para_k_otimo <- cumsum_k_monte_carlo %>%
  filter(k == 1 & desvio_detectavel == round(valor_otimo_k, 1)) %>%
  rename(`Φ₁` = h1_phi) %>%
  mutate(observacao = list(1:nH1))

ggplotly(
  cusum_para_k_otimo %>%
    unnest(
      cols = c(
        `Φ₁`,
        observacao,
        pos,
        neg,
        h1_residuos,
        fora_de_controle_pos,
        fora_de_controle_neg
      )
    ) %>%
    ggplot() +
    geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
    geom_hline(aes(
      yintercept = -5, color = "Limite de decisão"
    ), linetype = "dashed") +
    geom_line(aes(
      x = observacao, y = pos, color = "N+"
    )) +
    geom_line(aes(
      x = observacao, y = neg, color = "N-"
    )) +
    geom_point(
      data = . %>% filter(fora_de_controle_pos),
      aes(x = observacao, y = pos, color = "Fora de controle"),
      size = 1,
      shape = 4
    ) +
    geom_point(
      data = . %>% filter(fora_de_controle_neg),
      aes(x = observacao, y = neg, color = "Fora de controle"),
      size = 1,
      shape = 4
    ) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "N+" = adjustcolor("blue", alpha.f = 0.6),
        "N-" = adjustcolor("darkgreen", alpha.f = 0.6),
        "Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.6),
        "Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
      )
    ) +
    labs(title = "CUSUM para resíduos βAR(1)", ) +
    theme.base +
    theme(legend.position = "bottom", legend.title = element_blank()) +
    facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
  plotly.base

Resumo

Em média, diminuir o K aumenta a sensibilidade do CUSUM para detectar mudanças no processo.

cumsum_k_monte_carlo_resumo <- cumsum_k_monte_carlo %>%
  group_by(desvio_detectavel, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

datatable(
  cumsum_k_monte_carlo_resumo %>%
    mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e K, com H = 5",
  colnames = c("Valores de K", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)

Gráficos

Análise gráfica das médias de FPFCs.

ggplotly(
  cumsum_k_monte_carlo_resumo %>%
    ggplot(aes(
      x = h1_phi,
      y = mean,
      color = factor(desvio_detectavel)
    )) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      color = "Valores de K",
      title = "FPFCs por Φ₁ e K com H = 5"
    ) +
    theme.base
) %>%
  plotly.base

E, a distribuição dos PFCs.

ggplotly(
  cumsum_k_monte_carlo %>%
    ggplot(aes(
      x = factor(h1_phi),
      y = fracao_fora_de_controle,
      fill = factor(desvio_detectavel)
    )) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      fill = "Valores de K",
      title = "FPFCs por Φ₁ e K com H = 5"
    ) +
    theme.base
) %>%
  plotly.base %>%
  layout(boxmode = 'group')

Estimando o valor ótimo para H

Semelhante ao valor ótimo para K, queremos que, sob \(H_1\), nas 100 amostras iniciais a probabilidade de um ponto fora de controle seja de, nominalmente, 5%. Portanto, vamos buscar um valor para H onde isso acontece.

intervalo_de_decisao_otimo_optimize <- function(amostra_inicial) {
  for (x in seq(1, 6, by = 1)) {
    fora <- cusum_qcc(amostra_inicial, intervalo_de_decisao = x)$total_fora_de_controle
    if (fora == 5) {
      return(x)
    }
  }
  return(NA)
}

cusum_h_otimo <- cache_dados("cusum-h-otimo", \() {
  data.frame(id = 1:1000) %>%
    rowwise() %>%
    mutate(
      h0_amostra = list(barma.sim(
        n = nH0, phi = phi_parametro, seed = id
      )),
      h0_phi = list(barma.phi_estimado(h0_amostra)),
      h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
      minimo = intervalo_de_decisao_otimo_optimize(h0_residuos)
    )
})

ggplotly(
  cusum_h_otimo %>%
    filter(!is.na(minimo)) %>%
    ggplot(aes(
      x = 0,
      y = minimo,
      group = 0,
      fill = "red"
    )) +
    geom_boxplot() +
    labs(x = element_blank(), y = "Valor ótimo de K", title = "Valor ótimo de K por simulação") +
    theme.base +
    theme.no_legend
) %>%
  plotly.base

valor_otimo_h <- mean(cusum_h_otimo$minimo, na.rm = TRUE)

kable(
  cusum_h_otimo %>%
    filter(!is.na(minimo)) %>%
    group_by() %>%
    summarise(
      "Média" = mean(minimo),
      "Desvio padrão" = sd(minimo),
      "Amostras" = n(),
    ),
  caption = "Valor ótimo para H",
  align = 'c',
  digits = 3
)

Valor ótimo para H

Média	Desvio padrão	Amostras
3.16	0.712	125

Monte Carlo

cumsum_h_monte_carlo <- cache_dados("simulacao-cusum-h", function() {
  gerador_monte_carlo(parametros = list(intervalo_de_decisao = seq(1, 6, by = 1))) %>%
    mutate(
      qcc = list(
        cusum_qcc(
          h0_residuos,
          intervalo_de_decisao = intervalo_de_decisao,
          amostra_subsequente = h1_residuos
        )
      ),
      pos = list(qcc$pos),
      neg = list(qcc$neg),
      fora_de_controle_pos = list(qcc$fora_de_controle_pos),
      fora_de_controle_neg = list(qcc$fora_de_controle_neg),
      total_fora_de_controle = qcc$total_fora_de_controle,
      fracao_fora_de_controle = qcc$fracao_fora_de_controle
    ) %>%
    select(-qcc)
})

Cartas para K ótimo

Vamos analisar o CUSUM para o valor ótimo de K encontrado.

cusum_para_h_otimo <- cumsum_h_monte_carlo %>%
  filter(k == 1 & intervalo_de_decisao == round(valor_otimo_h)) %>%
  rename(`Φ₁` = h1_phi) %>%
  mutate(observacao = list(1:nH1))

ggplotly(
  cusum_para_h_otimo %>%
    unnest(
      cols = c(
        `Φ₁`,
        observacao,
        pos,
        neg,
        h1_residuos,
        fora_de_controle_pos,
        fora_de_controle_neg
      )
    ) %>%
    ggplot() +
    geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
    geom_hline(aes(
      yintercept = -5, color = "Limite de decisão"
    ), linetype = "dashed") +
    geom_line(aes(
      x = observacao, y = pos, color = "N+"
    )) +
    geom_line(aes(
      x = observacao, y = neg, color = "N-"
    )) +
    geom_point(
      data = . %>% filter(fora_de_controle_pos),
      aes(x = observacao, y = pos, color = "Fora de controle"),
      size = 1,
      shape = 4
    ) +
    geom_point(
      data = . %>% filter(fora_de_controle_neg),
      aes(x = observacao, y = neg, color = "Fora de controle"),
      size = 1,
      shape = 4
    ) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "N+" = adjustcolor("blue", alpha.f = 0.6),
        "N-" = adjustcolor("darkgreen", alpha.f = 0.6),
        "Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.8),
        "Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
      )
    ) +
    labs(title = "CUSUM para resíduos βAR(1)", ) +
    theme.base +
    theme(legend.position = "bottom", legend.title = element_blank()) +
    facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
  plotly.base

Resumo

Em média, diminuir o H aumenta a sensibilidade do CUSUM para detectar mudanças no processo.

cumsum_h_monte_carlo_resumo <- cumsum_h_monte_carlo %>%
  group_by(intervalo_de_decisao, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

datatable(
  cumsum_h_monte_carlo_resumo %>%
    mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e H, com K = 1",
  colnames = c("Valores de H", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)

Gráficos

Análise gráfica das médias de FPFCs.

ggplotly(
  cumsum_h_monte_carlo_resumo %>%
    ggplot(aes(
      x = h1_phi,
      y = mean,
      color = factor(intervalo_de_decisao)
    )) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      color = "Valores de H",
      title = "FPFCs por Φ₁ e H, com K = 1",
    ) +
    theme.base
) %>%
  plotly.base

E, a distribuição dos PFCs.

ggplotly(
  cumsum_h_monte_carlo %>%
    ggplot(aes(
      x = factor(h1_phi),
      y = fracao_fora_de_controle,
      fill = factor(intervalo_de_decisao)
    )) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      fill = "Valores de K",
      title = "FPFCs por Φ₁ e K com H = 5"
    ) +
    theme.base
) %>%
  plotly.base %>%
  layout(boxmode = 'group')

E se… juntarmos CUSUM(K) e CUSUM(H)?

Vamos verificar se existe alguma combinação ótima entre os valores de K e H para detectar pontos fora de controle.

comb_cusum <- cache_dados("simulacao-cusum-k-h", function() {
  gerador_monte_carlo(parametros = expand.grid(
    desvio_detectavel = seq(0.6, 1.8, by = 0.2),
    intervalo_de_decisao = seq(3, 6, by = 1)
  )) %>%
    mutate(
      qcc = list(
        cusum_qcc(
          h0_residuos,
          desvio_detectavel = desvio_detectavel,
          intervalo_de_decisao = intervalo_de_decisao,
          amostra_subsequente = h1_residuos
        )
      ),
      pos = list(qcc$pos),
      neg = list(qcc$neg),
      fora_de_controle_pos = list(qcc$fora_de_controle_pos),
      fora_de_controle_neg = list(qcc$fora_de_controle_neg),
      total_fora_de_controle = qcc$total_fora_de_controle,
      fracao_fora_de_controle = qcc$fracao_fora_de_controle
    ) %>%
    select(-qcc)
})

Resultados

comb_cusum_resumo <- comb_cusum %>%
  group_by(desvio_detectavel, intervalo_de_decisao, h1_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  ) %>%
  mutate(
    desvio_detectavel = as.factor(desvio_detectavel),
    intervalo_de_decisao = as.factor(intervalo_de_decisao),
  )

datatable(
  comb_cusum_resumo %>%
    mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle",
  colnames = c(
    "Valores de K",
    "Valores de H",
    "Valores de Φ₁",
    "Média",
    "Mínimo",
    "Máximo"
  )
)

ggplotly(
  comb_cusum_resumo %>%
    group_by(desvio_detectavel, intervalo_de_decisao) %>%
    summarise(
      h1_phi = list(h1_phi),
      mean = list(mean),
      .groups = "drop"
    ) %>%
    filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
    unnest(cols = c(h1_phi, mean)) %>%
    ggplot(
      aes(
        x = h1_phi,
        y = mean,
        color = desvio_detectavel,
        linetype = intervalo_de_decisao
      )
    ) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      color = "Valores de K",
      linetype = "Valores de H",
      title = "FPFCs por Φ₁, K e H"
    ) +
    theme.base
) %>%
  plotly.base

De uma forma geral, observa-se que aumentando K e H diminui a sensibilidade do CUSUM para detectar pontos fora de controle, e observa-se um leve aumento no poder quando combinamos H = 5 e K = 0.8, com um leve aumento na variabilidade.

EWMA x CUSUM(K/H)

Vamos comparar os algoritmos EWMA e CUSUM(K) e CUSUM(H) para diferentes valores de \(\Phi_1\).

Aqui chamamos de CUSUM(K) as estatísticas onde o \(H\) for mantido constante e o \(K\) é variado, e o inverso para CUSUM(H).

estatisticas_resumo_combinadas <- bind_rows(
  cumsum_k_monte_carlo_resumo %>%
    filter(desvio_detectavel == 0.9) %>%
    mutate(
      desvio_detectavel = as.factor(desvio_detectavel),
      algoritmo = "CUSUM(K)"
    ) %>%
    rename(parametro = desvio_detectavel),
  cumsum_h_monte_carlo_resumo %>%
    filter(intervalo_de_decisao == 4) %>%
    mutate(
      intervalo_de_decisao = as.factor(intervalo_de_decisao),
      algoritmo = "CUSUM(H)"
    ) %>%
    rename(parametro = intervalo_de_decisao),
  comb_cusum_resumo %>%
    filter(desvio_detectavel == 0.8 & intervalo_de_decisao == 5) %>%
    mutate(`(K,H)` = as.factor(
      paste(intervalo_de_decisao, desvio_detectavel, sep = ";")
    ), algoritmo = "CUSUM(K,H)") %>%
    rename(parametro = `(K,H)`),
  ewma_monte_carlo_resumo %>%
    filter(lambda == 0.7) %>%
    rename(parametro = lambda) %>% mutate(algoritmo = "EWMA(λ)")
) %>%
  select(algoritmo, parametro, h1_phi, mean) %>%
  mutate(parametro = as.factor(parametro))

ggplotly(
  estatisticas_resumo_combinadas %>%
    ggplot() +
    geom_line(aes(
      x = h1_phi,
      y = mean,
      color = parametro,
      linetype = algoritmo
    )) +
    geom_point(aes(
      x = h1_phi,
      y = mean,
      color = parametro,
      shape = algoritmo
    )) +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros",
      shape = "Algoritmo"
    ) +
    theme.base
) %>%
  plotly.base

estatisticas_combinadas <- bind_rows(
  cumsum_k_monte_carlo %>%
    filter(desvio_detectavel == 0.9) %>%
    mutate(
      desvio_detectavel = as.factor(desvio_detectavel),
      algoritmo = "CUSUM(K)"
    ) %>%
    rename(parametro = desvio_detectavel),
  cumsum_h_monte_carlo %>%
    filter(intervalo_de_decisao == 4) %>%
    mutate(
      intervalo_de_decisao = as.factor(intervalo_de_decisao),
      algoritmo = "CUSUM(H)"
    ) %>%
    rename(parametro = intervalo_de_decisao),
  comb_cusum %>%
    filter(desvio_detectavel == 0.8 & intervalo_de_decisao == 5) %>%
    mutate(`(K,H)` = as.factor(
      paste(intervalo_de_decisao, desvio_detectavel, sep = ";")
    ), algoritmo = "CUSUM(K,H)") %>%
    rename(parametro = `(K,H)`),
  ewma_monte_carlo %>%
    mutate(lambda = as.factor(lambda), algoritmo = "EWMA(λ)") %>%
    filter(lambda == 0.7) %>%
    rename(parametro = lambda)
) %>%
  select(algoritmo, parametro, h1_phi, fracao_fora_de_controle) %>%
  mutate(parametro = as.factor(parametro), h1_phi = as.factor(h1_phi))

ggplotly(
  estatisticas_combinadas %>%
    ggplot() +
    geom_boxplot(
      aes(
        x = h1_phi,
        y = fracao_fora_de_controle,
        fill = parametro,
        shape = algoritmo
      ),
    ) +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      fill = "Valores dos parâmetros",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros"
    ) +
    facet_wrap( ~ algoritmo) +
    theme.base
) %>%
  plotly.base %>%
  layout(boxmode = 'group')

Análise

A partir da análise da média dos pontos fora de controle, observa-se que o CUSUM(K), em média, possui poder semelhante ao CUSUM(H), enquanto o EWMA é inferior aos dois.

Já a análise gráfica dos boxplots mostra que o CUSUM(K) possui uma variabilidade, em geral, maior que o CUSUM(H), enquanto o EWMA possui, em geral, a menor variabilidade entre os algoritmos.

E se… combinarmos os dois algoritmos?

Vamos combinar EWMA e CUSUM e ver o que acontece.

comb_monte_carlo <- cache_dados("simulacao-ewma-cusum", function() {
  gerador_monte_carlo(
    parametros = expand.grid(
      desvio_detectavel = seq(0.2, 1.6, by = 0.2),
      intervalo_de_decisao = seq(1, 6, by = 1),
      lambda = seq(0.1, 1.0, by = 0.1)
    ),
    numero_de_execucoes = 50
  ) %>%
    mutate(
      "(K;H;λ)" = factor(
        paste(desvio_detectavel, intervalo_de_decisao, lambda, sep = ";")
      ),
      ewma = list(ewma_qcc(h0_residuos, lambda, h1_residuos)),
      cumsum = list(
        cusum_qcc(
          h0_residuos,
          desvio_detectavel = desvio_detectavel,
          intervalo_de_decisao = intervalo_de_decisao,
          amostra_subsequente = h1_residuos
        )
      ),
      fora_de_controle_coop = list(
        ewma$fora_de_controle |
          (cumsum$fora_de_controle_pos | cumsum$fora_de_controle_neg)
      ),
      fracao_fora_de_controle_coop = sum(fora_de_controle_coop) / nH1,
      fora_de_controle_opos = list(
        ewma$fora_de_controle &
          (cumsum$fora_de_controle_pos | cumsum$fora_de_controle_neg)
      ),
      fracao_fora_de_controle_opos = sum(fora_de_controle_opos) / nH1
    ) %>%
    select(-ewma, -cumsum)
})

EWMA + CUSUM

Vamos verificar como se comporta a combinação dos algoritmos EWMA e CUSUM, ou seja, se um dos algoritmos detectar um ponto fora de controle, então o processo é considerado fora de controle.

Nesta forma mantemos o nível de significância para ambos os algoritmos em 5%.

comb_coop_monte_carlo_resumo <- comb_monte_carlo %>%
  group_by(`(K;H;λ)`, h1_phi, desvio_detectavel, lambda) %>%
  summarise(
    mean = mean(fracao_fora_de_controle_coop),
    min = min(fracao_fora_de_controle_coop),
    max = max(fracao_fora_de_controle_coop),
    .groups = "drop"
  ) %>%
  group_by(`(K;H;λ)`) %>%
  summarise(
    h1_phi = list(h1_phi),
    desvio_detectavel = list(desvio_detectavel),
    lambda = list(lambda),
    mean = list(mean),
    min = list(min),
    max = list(max),
    .groups = "drop"
  ) %>%
  filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
  unnest(cols = c(h1_phi, desvio_detectavel, lambda, mean, min, max))

datatable(
  comb_coop_monte_carlo_resumo %>%
    select(-`(K;H;λ)`) %>%
    mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e (K;H;λ)",
  colnames = c("Valores de Φ₁", "K", "H", "λ", "Média", "Mínimo", "Máximo")
)

ggplotly(
  comb_coop_monte_carlo_resumo %>%
    group_by(`(K;H;λ)`) %>%
    summarise(
      h1_phi = list(h1_phi),
      mean = list(mean),
      .groups = "drop"
    ) %>%
    filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
    unnest(cols = c(h1_phi, mean)) %>%
    ggplot(aes(
      x = h1_phi, y = mean, color = `(K;H;λ)`
    )) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      color = "Constantes (K;H;λ)",
      title = "FPFC por Φ₁ e (K;H;λ) (coop)"
    ) +
    theme.base
) %>%
  plotly.base

EWMA x CUSUM

Vamos verificar como se comporta a combinação dos algoritmos EWMA vs. CUSUM, ou seja, apenas se ambos os algoritmos detectarem um ponto fora de controle, então o processo é considerado fora de controle.

Fazendo isso, diminui-se o nível de significância para cada um dos algoritmos. Ainda queremos que o nível de significância final seja de 5%, para tanto buscamos uma relação \(0.05 = 1 - (1 - \alpha_{EWMA})(1 - \alpha_{CUSUM})\), obtida a partir da significância para testes simultâneos.

comb_opos_monte_carlo_resumo <- comb_monte_carlo %>%
  group_by(`(K;H;λ)`, h1_phi, desvio_detectavel, lambda) %>%
  summarise(
    mean = mean(fracao_fora_de_controle_opos),
    min = min(fracao_fora_de_controle_opos),
    max = max(fracao_fora_de_controle_opos),
    .groups = "drop"
  ) %>%
  group_by(`(K;H;λ)`) %>%
  summarise(
    h1_phi = list(h1_phi),
    desvio_detectavel = list(desvio_detectavel),
    lambda = list(lambda),
    mean = list(mean),
    min = list(min),
    max = list(max),
    .groups = "drop"
  ) %>%
  filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
  unnest(cols = c(h1_phi, desvio_detectavel, lambda, mean, min, max))

datatable(
  comb_opos_monte_carlo_resumo %>%
    select(-`(K;H;λ)`) %>%
    mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
  caption = "Pontos fora de controle por Φ₁ e (K;H;λ)",
  colnames = c("Valores de Φ₁", "K", "H", "λ", "Média", "Mínimo", "Máximo")
)

ggplotly(
  comb_opos_monte_carlo_resumo %>%
    group_by(`(K;H;λ)`) %>%
    summarise(
      h1_phi = list(h1_phi),
      mean = list(mean),
      .groups = "drop"
    ) %>%
    filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.052)) %>%
    unnest(cols = c(h1_phi, mean)) %>%
    ggplot(aes(
      x = h1_phi, y = mean, color = `(K;H;λ)`
    )) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      color = "Constantes (K;H;λ)",
      title = "FPFC por Φ₁ e (K;H;λ) (opos)"
    ) +
    theme.base
) %>%
  plotly.base

A melhor combinação encontrada

estatisticas_combinadas_resumo_ew_sum <- bind_rows(
  estatisticas_resumo_combinadas,
  comb_coop_monte_carlo_resumo %>%
    filter(`(K;H;λ)` == "1.4;4;0.8") %>%
    mutate(parametro = as.factor(`(K;H;λ)`)) %>%
    mutate(algoritmo = "EW-SUM(coop)") %>%
    select(algoritmo, parametro, h1_phi, mean),
  comb_opos_monte_carlo_resumo %>%
    filter(`(K;H;λ)` == "1;4;0.1") %>%
    mutate(parametro = as.factor(`(K;H;λ)`)) %>%
    mutate(algoritmo = "EW-SUM(opos)") %>%
    select(algoritmo, parametro, h1_phi, mean)
)

ggplotly(
  estatisticas_combinadas_resumo_ew_sum %>%
    ggplot() +
    geom_line(aes(
      x = h1_phi,
      y = mean,
      color = parametro,
      linetype = algoritmo
    )) +
    geom_point(aes(
      x = h1_phi,
      y = mean,
      color = parametro,
      shape = algoritmo
    )) +
    labs(
      x = "Valores de Φ₁",
      y = "Fração de pontos fora de controle",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros",
      shape = "Algoritmo"
    ) +
    theme.base
) %>%
  plotly.base